////////////////////////////////////////////////////////////////////////
//irmm_Raster.h
//This is the description of classes of image library including
//======================================================================
//		Gaussian filtering
//		segmentation of color image using mean shift algorithm
//=======================================================================
//	Author: Zhuowen Tu
////////////////////////////////////////////////////////////////////////
#ifndef _IRMM_RASTER_H
#define _IRMM_RASTER_H

#include <assert.h>
#include <memory.h>
#include <math.h>
#include <stdlib.h>
#include <limits.h>
#include <fstream>
#include <strstream>
#include <memory.h>
#include <time.h>
#include <string.h>
#include <stdio.h>

#include "mcPublic.h"
#include "Matrix.h"
#include "Vector.h"
#include "mcLinkList.h"

#define		ORG_SIZE

#define Max_rects		50

#ifndef MIN
#define MIN(x,y)        ((x) <= (y) ? (x) : (y))
#endif

#define Min3(x,y,z)     (MIN((x),(y)) <= (z) ? MIN((x),(y)) : (z))

#define min(a, b)  (((a) < (b)) ? (a) : (b))
#define max(a, b)  (((a) > (b)) ? (a) : (b))

#define bzero(b,len)     memset(b,'\0',len)
#define bcopy(b1,b2,len) memmove(b2,b1,len)


#define nil 0


#define PPM 6
#define PGM 5
#define ASCIIPPM 3
#define ASCIIPGM 2
#define BINARYPPM 6
#define BINARYPGM 5


enum ImageFileFormat
{
	iffPpmAscii		= 3,
	iffPpmBinary	= 6,
	iffPgmAscii		= 2,
	iffPgmBinary	= 5
};


class irmm_Size
{
public:
	int		width;		// width of the window
	int		height;		// height of the window

public:
	// construction & destruction
	irmm_Size(){width=0; height=0;};
	irmm_Size(int h, int w){width=w; height=h;};

	// operators
	inline irmm_Size & operator=(const irmm_Size &a)
	{
		if (this != &a)
		{
			width = a.width;
			height = a.height;
		}

		return *this;
	}
	inline bool operator==(const irmm_Size &a)
	{
		return (height==a.height && width==a.width);
	}


	// operations
	void Set(int h, int w){width=w; height=h;};
};

// this is a general rectangle class that captures the bound
// of region
class irmm_Rect
{
public:
	int	left,right;
	int	top,bottom;

	// consturction & destruction
	irmm_Rect(){left=0;right=0;top=0;bottom=0;};
	irmm_Rect(const int l, const int t, const int r, const int b){left=l;right=r;top=t;bottom=b;};

	inline void SetRect(const int l, const int t, const int r, const int b){left=l;right=r;top=t;bottom=b;};
	//operators
	inline irmm_Rect & operator=(const irmm_Rect &a);
	inline bool operator==(const irmm_Rect &a);
	inline irmm_Rect operator&(const irmm_Rect a);
	inline irmm_Rect operator+(const McImagePoint a);

	inline  void	UnionRect(const irmm_Rect a, const irmm_Rect b);
	inline	bool	InRect(const int row, const int col) const;
	inline	bool	InRect(const irmm_Rect& rect) const;
	inline	bool	Overlap(const irmm_Rect& rect) const;
	inline	void	Set(int lLeft,int lTop,int lRight,int lBottom)
		{left=lLeft;top=lTop;right=lRight;bottom=lBottom;};
	inline	void	BoundRow(const int row);
	inline	void	BoundCol(const int col);
	inline	int	Width(void) const {return right-left+1;};
	inline	int	Height(void) const {return bottom-top+1;};
	inline	void	SetBound(int i, int j);
	inline	bool	ValidBound(void) const;
	inline	int	OverlappingArea(const irmm_Rect& rect) const;
	inline	int	Area(void) const;
	inline	void	Center(int &j, int &i) const;
	inline	void	NormalizeRect(void);
	inline	McImagePoint TopLeft(void);
};

irmm_Rect & irmm_Rect::operator=(const irmm_Rect &a)
{
	if (this != &a)
	{
		left	= a.left;
		right	= a.right;
		top		= a.top;
		bottom	= a.bottom;
	}

	return *this;
}

bool irmm_Rect::operator==(const irmm_Rect &a)
{
	return (left==a.left && right==a.right &&
			top==a.top && bottom==a.bottom);
}

inline irmm_Rect irmm_Rect::operator+(const McImagePoint a)
{
	irmm_Rect b;

	b.left = left+a.Col();
	b.top = top+a.Row();

	b.right = right+a.Col();
	b.bottom = bottom+a.Row();

	return b;
}


inline	McImagePoint irmm_Rect::TopLeft(void)
{
	McImagePoint a;
	a.SetRow(top);
	a.SetCol(left);

	return a;
}

inline void irmm_Rect::NormalizeRect(void)
{
	int temp;
	if (left > right)
	{
		temp = left;
		left = right;
		right = temp;
	}
	if (top > bottom)
	{
		temp = top;
		top = bottom;
		bottom = temp;
	}
}

inline void	irmm_Rect::SetBound(int i, int j)
{
	if (i<left)		left	= i;
	if (i>right)	right	= i;
	if (j<top)		top		= j;
	if (j>bottom)	bottom	= j;
}

inline bool	irmm_Rect::InRect(const int row, const int col) const
{
	if (row<top || row>bottom || col<left || col>right)
		return false;
	else
		return true;
}

inline	bool irmm_Rect::Overlap(const irmm_Rect& rect) const
{
	bool boverlap;
	boverlap = InRect(rect.top, rect.left) || InRect(rect.top, rect.right) ||
				InRect(rect.bottom, rect.left) || InRect(rect.bottom, rect.right);
	if (boverlap)
		return true;
	else
	{
		boverlap =  (rect.InRect(top, left) || rect.InRect(top, right) ||
					rect.InRect(bottom, left) || rect.InRect(bottom, right));
		if (boverlap)
			return true;
		else
		{
			McPoint2D ppts1[5],ppts2[5];
			ppts1[0].x = top;		ppts1[0].y = left;
			ppts1[1].x = top;		ppts1[1].y = right;
			ppts1[2].x = bottom;	ppts1[2].y = right;
			ppts1[3].x = bottom;	ppts1[3].y = left;
			ppts1[4] = ppts1[0];

			ppts2[0].x = rect.top;		ppts2[0].y = rect.left;
			ppts2[1].x = rect.top;		ppts2[1].y = rect.right;
			ppts2[2].x = rect.bottom;	ppts2[2].y = rect.right;
			ppts2[3].x = rect.bottom;	ppts2[3].y = rect.left;
			ppts2[4] = ppts2[0];

			for (int i=0; i<4; i++)
				for (int j=0; j<4; j++)
					if (LineIntersection(NULL,ppts1[i],ppts1[i+1],ppts2[j],ppts2[j+1])==CROSSING)
						return true;
			return false;
		}
	}
}

inline bool	irmm_Rect::InRect(const irmm_Rect& rect) const
{
	int row1=rect.top,		col1=rect.left;
	int row2=rect.bottom,	col2=rect.left;
	int row3=rect.top,		col3=rect.right;
	int row4=rect.bottom,	col4=rect.right;

	if (InRect(row1,col1) && InRect(row2,col2) && InRect(row3,col3) && InRect(row4,col4))
		return true;
	else
		return false;
}

inline	void irmm_Rect::BoundRow(const int row)
{
	if (row<top) top = row;
	if (row>bottom) bottom = row;
}

inline	void irmm_Rect::BoundCol(const int col)
{
	if (col<left) left = col;
	if (col>right) right = col;
}

inline bool irmm_Rect::ValidBound(void) const
{
	if (left<=right && top<=bottom)
		return true;
	else
		return false;
}

inline int irmm_Rect::Area(void) const
{
	if (ValidBound())
		return ((int)Height())*Width();
	else
		return 0;
}

inline	void irmm_Rect::Center(int &j, int &i) const
{
	j = (top+bottom)/2;
	i = (left+right)/2;
}

inline int irmm_Rect::OverlappingArea(const irmm_Rect& rect) const
{
	int row_top,row_bottom,col_left,col_right;

	row_top		= max(top, rect.top);
	row_bottom	= min(bottom, rect.bottom);
	col_left	= max(left, rect.left);
	col_right	= min(right, rect.right);

	if (row_top<=row_bottom && col_left<=col_right)
		return (((int)(row_bottom-row_top+1))*(col_right-col_left+1));
	else
		return 0;
}

inline  void irmm_Rect::UnionRect(const irmm_Rect a, const irmm_Rect b)
{
	top		= min(b.top, a.top);
	bottom	= max(b.bottom, a.bottom);
	left	= min(b.left, a.left);
	right	= max(b.right, a.right);
}


inline irmm_Rect irmm_Rect::operator&(const irmm_Rect a)
{
	irmm_Rect o;

	o.top		= max(top, a.top);
	o.bottom	= min(bottom, a.bottom);
	o.left		= max(left, a.left);
	o.right		= min(right, a.right);

	if (o.top > o.bottom)
		o.bottom = o.top-1;

	if (o.left > o.right)
		o.right = o.left-1;

	return o;
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
class MCMC_PixelSet
{
public:
	MCMC_PixelSet(){;};
	~MCMC_PixelSet(){;};
	
	inline MCMC_PixelSet & operator=(MCMC_PixelSet &a);

	inline	void Empty(void){list_pixels_row.Empty(); list_pixels_col.Empty();};
	inline	void Insert(const int &j, const int &i){list_pixels_row.Insert(j); list_pixels_col.Insert(i);};
	inline	void InsertToHead(const int &j, const int &i){list_pixels_row.InsertToHead(j); list_pixels_col.InsertToHead(i);};
	inline	bool	GetFirst(int &j, int &i){list_pixels_row.GetFirst(j); return list_pixels_col.GetFirst(i);};
	inline	bool	GetNext(int &j, int &i){list_pixels_row.GetNext(j); return list_pixels_col.GetNext(i);};
	inline	bool	DeleteTheFirstPixel(void){list_pixels_row.DeleteTheFirst(); return list_pixels_col.DeleteTheFirst();};
	inline	int		Num(void){return list_pixels_row.Num();};
	inline	void	AddFrom(MCMC_PixelSet &pixel_set);
	
	inline	bool	GetLast(int &j, int &i){list_pixels_row.GetLast(j); return list_pixels_col.GetLast(i);};
	inline	bool	GetPrev(int &j, int &i){list_pixels_row.GetPrev(j); return list_pixels_col.GetPrev(i);};
	inline	bool	DeleteTheLastPixel(void){list_pixels_row.DeleteTheLast(); return list_pixels_col.DeleteTheLast();};
	
	inline	irmm_Rect	Bound(void);

public:
	double	energy;

protected:
	LinkList<int>	list_pixels_row;
	LinkList<int>	list_pixels_col;
};

typedef MCMC_PixelSet* PMCMC_PixelSet;

inline void	MCMC_PixelSet::AddFrom(MCMC_PixelSet &pixel_set)
{
	bool	bgo;
	int		i,j;

	bgo = pixel_set.GetFirst(j, i);
	while (bgo)
	{
		Insert(j, i);
		bgo = pixel_set.GetNext(j, i);
	}
}

inline MCMC_PixelSet & MCMC_PixelSet::operator=(MCMC_PixelSet &a)
{
	if (this != &a)
	{
		list_pixels_row = a.list_pixels_row;
		list_pixels_col = a.list_pixels_col;

		energy = a.energy;
	}

	return *this;
}

inline irmm_Rect MCMC_PixelSet::Bound(void)
{
	bool	bgo;
	int		i,j;
	irmm_Rect	rect;

	bgo = GetFirst(j, i);
	if (bgo)
		rect.Set(i, j, i, j);
	while (bgo)
	{
		rect.SetBound(i, j);
		
		bgo = GetNext(j, i);
	}

	return rect;
}


///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////


#define GAUSS 0
#define FIRST_DERIV 1
#define SECOND_DERIV 2

static int row_neighbors[8]={1,0,-1,-1,-1,0,1,1};
static int col_neighbors[8]={-1,-1,-1,0,1,1,1,0};

static int row_neighbors4[4]={0,-1,0,1};
static int col_neighbors4[4]={-1,0,1,0};

static int row_neighbors4n[4]={1,-1,-1,1};
static int col_neighbors4n[4]={-1,-1,1,1};

/////////////////////////////////////////////////////////
// this is the class of raster image in one channel
template<class Tp>
class Raster: public mcmcMatrix<Tp>
{
public:
    // construction
	Raster();
	Raster(const Raster<Tp>& r);
	Raster(const int rows, const int cols);
    ~Raster();

	// operators
	Raster<Tp>&		operator=(const Raster<Tp>& r);
	Raster&			operator=(const mcmcMatrix<Tp>& r);
	Raster&			operator|(const Raster& r);
	Raster&			operator&(const Raster& r);

	// access
    Tp				Get(const int Row, const int Col) const;
	void            Set(Tp value);
	void			Set(const int Row, const int Col, Tp value);
	inline int		Width() const {return cols();};
	inline int		Height()const {return rows();};
	inline void		Zero(void){InitValue(0);};

	// operations
	virtual Raster<uchar> ToVisibleRaster(const bool bBinary=false);

	// merge dangling pixels
	virtual	void	MergeDanglingPixels(void);
	// used after diffusion
	virtual int		Relabel(const bool bMerge=true, mcmcVector<int> *pvc_original_id=NULL);
	// used for cleaning proposals for split
	virtual int		RelabelAndMergeTinyRegions(const int min_size=3);
	// used for an initial label map (including 0)
	virtual int		RelabelReorder(const int label_start_from=0);

	virtual bool	DanglingPixel(const int row, const int col);
	virtual int		FindUnion(mcmcMatrix<int> &mx_representative, mcmcMatrix<int> &mx_num_pixels);
	virtual int		GetRepresentative(mcmcMatrix<int> &mx_representative, const int j, const int i,int &m, int &n);

	virtual void	CopySubImage(irmm_Rect rect1, Raster<Tp>& image_dest, irmm_Rect rect2);
	virtual void	CopySubImage(irmm_Rect rect1, Raster<Tp>& image_dest);
	virtual void	ResizeTo(Raster<Tp>& ras_new_image, const char flag);

	virtual void	DrawRect(const irmm_Rect &rc, const Tp intensity);

private:
	inline  double	R(const double x);
	virtual void	DrawRectA(const irmm_Rect& rect, const Tp intensity);
};

extern const double     Un_prime;
extern const double     Vn_prime;
extern const double     Lt;

// Constants for LUV transformation
extern const double     Xn;
extern const double     Yn;
extern const double     Zn;

template<class Tp1, class Tp2> void CopyRaster(Raster<Tp1>& ras_dest, const Raster<Tp2>& ras_src);



//==========================================================
// Class : Raster
//==========================================================

//-----------------------------------------------------------
// construction and destruction

/////////////////////////////////////////////////////////////
template<class Tp>
Raster<Tp>::Raster() // Jiadi
{
}

template<class Tp>
Raster<Tp>::Raster(const Raster<Tp>& r) // Jiadi
{
	(*this) = r;
}

template<class Tp>
Raster<Tp>::Raster(const int rows, const int cols) // Jiadi
{
	Raster();

	SetDimension(rows, cols);
}


/////////////////////////////////////////////////////////////
template<class Tp>
Raster<Tp>::~Raster() // Jiadi
{
}

//-----------------------------------------------------------
// operators
/////////////////////////////////////////////////////////////
template<class Tp>
Raster<Tp>&	Raster<Tp>::operator=(const Raster<Tp>& r) // Jiadi
{
	*((mcmcMatrix<Tp>*) this) = *((mcmcMatrix<Tp>*) &r);

	return *this;
}

template<class Tp>
Raster<Tp>&	Raster<Tp>::operator=(const mcmcMatrix<Tp>& r) // Jiadi
{
	*((mcmcMatrix<Tp>*) this) = *((mcmcMatrix<Tp>*) &r);

	return *this;
}

template<class Tp1, class Tp2>
void CopyRaster(Raster<Tp1>& ras_dest,const Raster<Tp2>& ras_src) // Jiadi
{
	CopyMatrix(*((mcmcMatrix<Tp1>*) &ras_dest),*((mcmcMatrix<Tp2>*) &ras_src));
}

////////////////////////////////////////////////////////////////
template<class Tp>
Raster<Tp>& Raster<Tp>::operator|(const Raster<Tp>& r) // Jiadi ! questions about type casting <uchar>
{
	if (this != &r && rows() == r.rows() && cols() == r.cols())
		bitwise_or(data, r.data, data);
	return *this;
}

////////////////////////////////////////////////////////////////
template<class Tp>
Raster<Tp>& Raster<Tp>::operator&(const Raster<Tp>& r) // Jiadi
{
	if (this != &r && rows() == r.rows() && cols() == r.cols())
		bitwise_and(data, r.data, data);
	return *this;
}


//-----------------------------------------------------------
// access

////////////////////////////////////////////////////////////////
template<class Tp>
Tp Raster<Tp>::Get(const int Row, const int Col) const // Jiadi
{
	return Data(Row, Col);
}

////////////////////////////////////////////////////////////////
template<class Tp>
void Raster<Tp>::Set(const int Row, const int Col, Tp value) // Jiadi
{
	if (Row < rows() && Col < cols())
		(*this)(Row, Col) = value;
}

template<class Tp>
void Raster<Tp>::Set(Tp value)						 // Jiadi
{
	for (int i = 0; i < rows();  i++)
		for (int j = 0; j < cols();  j++)
			(*this)(i,j) = value;
}

template<class Tp>
void Raster<Tp>::CopySubImage(irmm_Rect rect1, Raster<Tp>& image_dest, irmm_Rect rect2) // Jiadi
{
	int i,j,l,m;

	m = rect2.top-rect1.top;
	l = rect2.left-rect1.left;

	for (j=rect1.top; j<=rect1.bottom; j++)
		for (i=rect1.left; i<=rect1.right; i++)
			image_dest(j+m, i+l) = Data(j,i);
}

template<class Tp>
void Raster<Tp>::CopySubImage(irmm_Rect rect1, Raster<Tp>& image_dest) // Jiadi
{
	int i,j;

	image_dest.SetDimension(rect1.Height(), rect1.Width());

	for (j=rect1.top; j<=rect1.bottom; j++)
		for (i=rect1.left; i<=rect1.right; i++)
			image_dest(j-rect1.top, i-rect1.left) = Data(j,i);
}

// transform current image to general image in char that could be seen
// directly 0~255
template<class Tp>
Raster<uchar> Raster<Tp>::ToVisibleRaster(const bool bBinary) // Jiadi
{
	Raster<uchar>	I;
	Tp				v1,v2;
	int			i;

	I.SetDimension(rows(), cols());
	GetMinMax(v1, v2);
	if (v1 != v2)
	{
		if (!bBinary)
			for (i=0; i<size(); i++)
				I(i) = uchar((double(Data(i)-v1))/(v2-v1)*255);
		else
			// as binary image
		{
			for (i=0; i<size(); i++)
			{
				if (Data(i)>0)
					I(i) = 255;
				else
					I(i) = 0;
			}
		}
	}
	return I;
}

inline bool IsDigit(const int a) // Jiadi
{
	if (a>='0' && a<='9')
		return true;
	else
		return false;
}


template<class Tp>
bool Raster<Tp>::DanglingPixel(const int row, const int col) // Jiadi
{
	int		i, m, n;
	bool	bDangling = false;

	if (Valid(row, col))
	{
		bDangling = true;
		for (i = 0; i<4; i++)
		{
			n = row + row_neighbors4[i];
			m = col + col_neighbors4[i];
			if (Valid(n, m) && Data(row, col) == Data(n, m))
			{
				bDangling = false;
				break;
			}
		}
		int temp = 0;
		if (bDangling)
		{
			for (i = 0; i<4; i++)
			{
				n = row + row_neighbors4n[i];
				m = col + col_neighbors4n[i];
				if (Valid(n, m) && Data(row, col) == Data(n, m))
					temp++;
			}
			if (temp<2)
				bDangling = false;
		}
	}
	return bDangling;
}


template<class Tp>
void Raster<Tp>::MergeDanglingPixels(void) // Jiadi
{
	int			j, i, m, n;

	for (j = 0; j<rows(); j++)
	for (i = 0; i<cols(); i++)
	{
		m = j + 1; n = i - 1;
		if (Valid(m, n))
		{
			if (Data(m, n) == Data(j, i) && Data(m, i) != Data(j, i) && Data(j, n) != Data(j, i))
			{
				if (!DanglingPixel(m, i))
				{
					(*this)(m, i) = Data(j, i);
					assert(Data(m, i) != Data(j, n));
				}
			}
		}
		m = j + 1; n = i + 1;
		if (Valid(m, n))
		{
			if (Data(m, n) == Data(j, i) && Data(m, i) != Data(j, i) && Data(j, n) != Data(j, i))
			{
				if (!DanglingPixel(m, i))
				{
					(*this)(m, i) = Data(j, i);
					assert(Data(m, i) != Data(j, n));
				}
			}
		}
	}
}

template<class Tp>
int Raster<Tp>::GetRepresentative(mcmcMatrix<int> &mx_representative, const int j, const int i,
								   int &m, int &n)  // Jiadi
{
	int k,l,temp;

	k = j; l = i;
	while (mx_representative.Data(k,l)>=0)
	{
		// get its parent
		temp = mx_representative.Data(k,l);
		// get the corresponding position
		mx_representative.GetPos(temp,k,l);
	}
	m = k; n = l;
	return mx_representative.Data(m,n);
}

// return: the largest value in the matrix
//  mx_representative: the representative for each pixel
//  mx_num_pixels: the total number pixels for each group
template<class Tp>
int Raster<Tp>::FindUnion(mcmcMatrix<int> &mx_representative, mcmcMatrix<int> &mx_num_pixels)  // Jiadi
{
	int	j,i,j1,i1,j2,i2,m,n,m1,n1,m2,n2;
	int		largest_label=-1;

	// initialization
	mx_representative.SetDimension(rows(),cols());
	mx_representative.InitValue(-1);
	mx_num_pixels.SetDimension(rows(),cols());
	mx_num_pixels.InitValue(1);

	for (j=0; j<rows(); j++)
	{
		for (i=0; i<cols(); i++)
		{
			if ((Tp)largest_label<Data(j,i))
				largest_label=(int)Data(j,i);

			m = j; n = i;
			// upper neighbor
			j1 = j-1;	i1 = i;
			if (Valid(j1,i1) && Data(j1,i1)==Data(j,i))
			{
				GetRepresentative(mx_representative,j1,i1,m1,n1);
				mx_representative(j,i) = m1*cols()+n1;
				mx_representative(m1,n1) = -max(abs(mx_representative(m1,n1)), 2);
				mx_num_pixels(j,i) = 0;
				mx_num_pixels(m1,n1) = mx_num_pixels(m1,n1)+1;

				m = m1; n = n1;
			}
			// left neighbor
			j2 = j;		i2 = i-1;
			if (Valid(j2,i2) && Data(j2,i2)==Data(j,i))
			{
				GetRepresentative(mx_representative,j2,i2,m2,n2);
				if (!(m2==m && n2==n))
				{
					if (abs(mx_representative(m,n))>abs(mx_representative(m2,n2)))
					{
						mx_representative(m,n) = -max(abs(mx_representative(m,n)),abs(mx_representative(m2,n2))+1);
						mx_representative(m2,n2) = m*cols()+n;
						mx_num_pixels(m,n) = mx_num_pixels(m,n)+mx_num_pixels(m2,n2);
						mx_num_pixels(m2,n2) = 0;
					}
					else
					{
						mx_representative(m2,n2) = -max(abs(mx_representative(m2,n2)),abs(mx_representative(m,n))+1);
						mx_representative(m,n) = m2*cols()+n2;
						mx_num_pixels(m2,n2) = mx_num_pixels(m,n)+mx_num_pixels(m2,n2);
						mx_num_pixels(m,n) = 0;
					}
				}
			}
		}
	}
	return largest_label;
}


#define MAX_REGION_NUM  2000 //1200

#define MAX_REGION_SIZENOTMERGED 1000 //900

#define MERGE_DANGLING_PIXELS				true
#define DONOTMERGE_DANGLING_PIXELS			false
#define MERGE_DISCONNECTED_REGIONS			true
#define DONOTMERGE_DISCONNECTED_REGIONS		false

template<class Tp>
int	Raster<Tp>::Relabel(const bool bMerge, mcmcVector<int> *pvc_original_id) //Jiadi
{
	mcmcMatrix<int>	mx_representative,mx_num_pixels;
	int			j,i,m,n,m1,n1,index;
	mcmcMatrix<int>		delsign;
	int				k,largest_label;
	mcmcVector<int>		rep_row(MAX_REGION_NUM),rep_col(MAX_REGION_NUM),num(MAX_REGION_NUM);
	bool			bOK=false;
	int				size_current_region;

	// initialization
	MergeDanglingPixels();
	if (pvc_original_id != NULL)
	{
		pvc_original_id->SetRowVector(MAX_REGION_NUM);
		pvc_original_id->InitValue(-1);
	}

	while (!bOK)
	{
		bOK = true;

		// mx_representative keeps the position of parent of pixel at (j,i)
		// mx_num_pixels keeps the number of connected pixles with the same id
		largest_label = FindUnion(mx_representative,mx_num_pixels);

		delsign.SetDimension(rows(),cols());
		delsign.InitValue(-1);

		for (i=0; i<MAX_REGION_NUM; i++)
			num(i) = 0;

		// find which region should be kept and which one should be merged
		for (j=0; j<rows(); j++)
			for (i=0; i<cols(); i++)
			{
				if (mx_representative(j,i)<0)
				{
					// get region id
					k = (int)Data(j,i);
					if (k>=MAX_REGION_NUM)
						return MAX_REGION_NUM+1;

					if (abs(mx_num_pixels(j,i))>num(k))
					{
						// keep this region
						rep_row(k) = j;
						rep_col(k) = i;
						num(k) = abs(mx_num_pixels(j,i));
					}
				}
			}
		for (k=0; k<MAX_REGION_NUM; k++)
		{
			if (num(k)>0)
				delsign(rep_row(k),rep_col(k))=-2;
		}

		for (j=0; j<rows(); j++)
			for (i=0; i<cols(); i++)
			{
				// find its representative
				GetRepresentative(mx_representative,j,i,m,n);
				index = (int)Data(m,n);

				if (delsign(m,n)==-1)
				{
					bOK = false;
					size_current_region = mx_num_pixels(m,n);

					if (bMerge && size_current_region<MAX_REGION_SIZENOTMERGED)
					{
						bool bfound=false;
						// find a neighbor to which this region merges
						for (k=0; k<4; k++)
						{
							m1 = m+row_neighbors4[k];
							n1 = n+col_neighbors4[k];
							if (Valid(m1,n1) && Data(m1,n1)!=Data(m,n))
							{
								int m2,n2;
								GetRepresentative(mx_representative,m1,n1,m2,n2);
								if (delsign(m2,n2)==-2)
								{
									bfound = true;
									break;
								}
							}
						}
						if (!bfound)
						{
							for (int jj=0; jj<rows(); jj++)
							{
								for (int ii=0; ii<cols(); ii++)
								{
									int mm1,nn1;
									GetRepresentative(mx_representative,jj,ii,mm1,nn1);
									if (mm1==m && nn1==n)
									{
										for (k=0; k<4; k++)
										{
											m1 = jj+row_neighbors4[k];
											n1 = ii+col_neighbors4[k];
											if (Valid(m1,n1) && Data(m1,n1)!=Data(m,n))
											{
												int m2,n2;
												GetRepresentative(mx_representative,m1,n1,m2,n2);
												if (delsign(m2,n2)==-2)
												{
													bfound = true;
													break;
												}
											}
										}
									}
									if (bfound) break;
								}
								if (bfound) break;
							}
						}
						if (bfound)
							delsign(m,n) = (int)Data(m1,n1);
					}
					else
					{
						// assign this region with a new id
						delsign(m,n) = (int)(++largest_label);

						// save its original id
						if (pvc_original_id != NULL)
							(*pvc_original_id)(delsign(m,n)) = (int)((*this)(m,n));
					}
					if (delsign(m,n)>=0)
						(*this)(j,i) = (Tp)delsign(m,n);
				}
				else if (delsign(m,n)>=0)
					(*this)(j,i) = (Tp)delsign(m,n);
			}
	}
	return largest_label;
}


template<class Tp>
int	Raster<Tp>::RelabelAndMergeTinyRegions(const int min_size)  //Jiadi
{
	mcmcMatrix<int>	union1,union2;
	int				num_regions,k;
	int			num_total_pixel,i,j,m,n;
	mcmcMatrix<int>		mx_index;

	mx_index.SetDimension(rows(), cols());
	mx_index.InitValue(0);

	num_total_pixel = 0;
	num_regions = 0;
	// find connected regions
	FindUnion(union1, union2);
	for (j=0; j<rows(); j++)
		for (i=0; i<cols(); i++)
		{
			if (Data(j,i)>0)
			{
				GetRepresentative(union1,j,i,m,n);
				if (mx_index(m, n)==0)
					mx_index(m, n) = ++num_regions;
				(*this)(j, i) = (Tp)mx_index(m, n);

				num_total_pixel++;
			}
		}

	if (num_regions <= 1)
		return num_regions;
	if (num_total_pixel <= min_size)
	{
		// assign each pixel with the value of 1
		for (j=0; j<rows(); j++)
			for (i=0; i<cols(); i++)
				if (Data(j,i)>0)
					(*this)(j, i) = 1;
		return 1;
	}


	mcmcMatrix<int>	mx_nb_topology,mx_rlabel;
	int			lb1,lb_nbs[2];
	mcmcVector<int> vc_counter,vc_label_final;

	vc_counter.SetRowVector(num_regions+1);
	vc_counter.InitValue(0);

	mx_nb_topology.SetDimension(num_regions+1, num_regions+1);
	mx_nb_topology.InitValue(0);
	vc_label_final.SetDimension(num_regions+1, 1);

	// initialize all the bounds
	for (i=0; i<vc_label_final.size(); i++)
		vc_label_final(i) = i;

	// count numbers and set relationships
	for (j=0; j<rows(); j++)
		for (i=0; i<cols(); i++)
		{
			if (Data(j,i)>0)
			{
				lb1 = (int)Data(j, i);
				// update region features
				lb_nbs[0] = -1;
				// check its right side
				if (Valid(j,i+1) && Data(j, i+1)>0)
					lb_nbs[0] = (int)Data(j, i+1);
				lb_nbs[1] = -1;
				// check the pixel underneath
				if (Valid(j+1, i) && Data(j+1, i)>0)
					lb_nbs[1] = (int)Data(j+1,i);
				for (k=0; k<2; k++)
				{
					// mark as neighbors
					if (lb_nbs[k]>0 && lb_nbs[k]!=lb1)
					{
						mx_nb_topology(lb1,lb_nbs[k]) = 1;
						mx_nb_topology(lb_nbs[k],lb1) = 1;
					}
				}
				vc_counter(lb1)++;
			}
		}

	bool	bvalid_vall=false;
	int		j1;

	j1 = 0;
	while (!bvalid_vall)
	{
		bvalid_vall = true;
		for (j=j1; j<vc_counter.size(); j++)
		{
			if (vc_counter(j)>0 && vc_counter(j)<min_size)
			{
				bvalid_vall = false;
				break;
			}
			else
				j1 = j;
		}
		if (!bvalid_vall)
		{
			for (i=0; i<vc_counter.size(); i++)
			{
				if (mx_nb_topology(j,i)>0 && vc_counter(i)>0 && i!=j)
				{
					vc_counter(i) += vc_counter(j);
					vc_counter(j) = 0;
					for (k=0; k<vc_counter.size(); k++)
					{
						if (vc_label_final(k)==j)
							vc_label_final(k) = i;
						if (mx_nb_topology(j,k)>0)
						{
							mx_nb_topology(i,k) = 1;
							mx_nb_topology(k,i) = 1;
						}
					}

					num_regions--;
					break;
				}
			}
			if (i >= vc_counter.size())
				vc_counter(j) = min_size;
		}
	}
	// rewrite labels
	for (j=0; j<rows(); j++)
		for (i=0; i<cols(); i++)
		{
			lb1 = (int)Data(j,i);
			if (lb1>0 && vc_counter(lb1)==0)
				// update region features
				(*this)(j,i) = (Tp)vc_label_final(lb1);
		}

	num_regions = 0;
	for (i=0; i<vc_counter.size(); i++)
	{
		if (vc_counter(i)>0)
		{
			num_regions++;
			vc_label_final(i) = num_regions;
		}
	}

	// change the labels
	for (j=0; j<rows(); j++)
		for (i=0; i<cols(); i++)
		{
			lb1 = (int)Data(j,i);
			if (lb1 > 0)
				// update region features
				(*this)(j,i) = (Tp)vc_label_final(lb1);
		}
	return num_regions;
}


template<class Tp>
int	Raster<Tp>::RelabelReorder(const int label_start_from) // Jiadi
{
	int *pindex,index,i,j,nmax,nmin,npossible,temp;
	
	nmax = -99999;
	nmin = 99999;
	
	for (j=0; j<rows(); j++)
		for (i=0; i<cols(); i++)
		{
			temp = (int)Data(j, i);
			if (temp >= label_start_from)
			{
				if (temp < nmin)	nmin = temp;
				if (temp > nmax)	nmax = temp;
			}
		}

	npossible = nmax+1-nmin;
	if (npossible <= 0)
		return 0;

	pindex = new int[npossible];
	for (i=0; i<npossible; i++)
		pindex[i] = -1;

	for (j=0; j<rows(); j++)
		for (i=0; i<cols(); i++)
		{
			temp = (int)Data(j, i);
			if (temp >= label_start_from)
				pindex[(temp-nmin)] = 1;
		}

	index = 1;
	for (i=0; i<npossible; i++)
	{
		if (pindex[i]>0)
		{
			pindex[i] = index;
			index++;
		}
	}
	for (j=0; j<rows(); j++)
		for (i=0; i<cols(); i++)
		{
			temp = (int)Data(j, i);
			if (temp >= label_start_from)
				(*this)(j, i) = (Tp)pindex[temp-nmin];
		}

	delete []pindex;

	return index-1;
}

#define TruncateNegative(a) (a >= 0.0 ? a : 0.0)

// flag:
#define INTRP_BILINEAR 0
#define INTRP_BICUBIC  1

template<class Tp>
inline double Raster<Tp>::R(const double x) // Jiadi
{
	double y,d=0.0;

	y = TruncateNegative(x+2);
	d = y*y*y;

	y = TruncateNegative(x+1);
	d += -4.0*y*y*y;

	y = TruncateNegative(x);
	d += 6.0*y*y*y;

	y = TruncateNegative(x-1);
	d += -4.0*y*y*y;

	return d/6.0;
}

template<class Tp>
void Raster<Tp>::ResizeTo(Raster<Tp>& ras_new, const char flag) // Jiadi
{
	double	w1,h1,w2,h2,x,y,dx,dy;
	int		j1,i1,j2,i2,row,col,m,n;
	Tp		f00,f01,f10,f11;

	if (size()<=1 || ras_new.size()<=1)
		return;

	// initialization
	h1 = rows()-1;			w1 = cols()-1;
	h2 = ras_new.rows()-1;	w2 = ras_new.cols()-1;
	h2 = h1/h2;				w2 = w1/w2;

	switch (flag)
	{
	case INTRP_BILINEAR:

			for (row=0; row<ras_new.rows(); row++)
				for (col=0; col<ras_new.cols(); col++)
				{
					x = h2*row;
					y = w2*col;

					j1 = (int)(floor(x));	i1 = (int)(floor(y));
					j2 = (int)(ceil(x));	i2 = (int)(ceil(y));

					j2 = j2 <= nrow-1 ? j2 : nrow-1;
					i2 = i2 <= ncol-1 ? i2 : ncol-1;

					dx = x-j1;	dy = y-i1;

					f00 = Data(j1, i1); f01 = Data(j1, i2);
					f10 = Data(j2, i1); f11 = Data(j2, i2);

					ras_new(row, col) = (Tp)((f10-f00)*dx + (f01-f00)*dy + (f11+f00-f01-f10)*dx*dy + f00);
				}
			break;

	case INTRP_BICUBIC:

			double temp;

			for (row=0; row<ras_new.rows(); row++)
				for (col=0; col<ras_new.cols(); col++)
				{
					x = h2*row;
					y = w2*col;

					j1 = (int)floor(x);	i1 = (int)floor(y);

					dx = x-j1;	dy = y-i1;

					temp = 0;
					for (m=-1; m<=2; m++)
						for (n=-1; n<=2; n++)
							if (Valid(j1+m, i1+n))
								temp += R(-dx+m) * R(dy-n) * Data(j1+m, i1+n);

					ras_new(row, col) = (Tp)temp;
				}
			break;
	}
}

template<class Tp>
void Raster<Tp>::DrawRectA(const irmm_Rect& rect, const Tp intensity) // Jiadi
{
	int			x,y;
	int			w = Width();
	int			h = Height();
	irmm_Rect	rcImg;

	rcImg.left = 0;	rcImg.top = 0;
	rcImg.bottom = rows()-1; rcImg.right = cols();

	for (x=rect.left;x<=rect.right;x++)
	{
		if (rcImg.top <= rect.top && x < rcImg.right)
			(*this)(rect.top, x) = intensity;
		if (rcImg.bottom > rect.bottom && x < rcImg.right)
			(*this)(rect.bottom, x) = intensity;
	}
	for (y=rect.top+1;y<rect.bottom;y++)
	{
		if (rcImg.left <= rect.left && y < rcImg.bottom)
			(*this)(y, rect.left) = intensity;
		if (rcImg.right > rect.right && y < rcImg.bottom)
			(*this)(y, rect.right) = intensity;
	}
}

template<class Tp>
void Raster<Tp>::DrawRect(const irmm_Rect& rc, const Tp intensity) // Jiadi
{
	irmm_Rect rect(rc.left-1, rc.top-1, rc.right+1, rc.bottom+1);
	DrawRectA(rect, intensity);
	rect = irmm_Rect(rc.left-2, rc.top-2, rc.right+2, rc.bottom+2);
	DrawRectA(rect, intensity);
	rect = irmm_Rect(rc.left-3, rc.top-3, rc.right+3, rc.bottom+3);
	DrawRectA(rect, intensity);
}

///////////////////////////   Supporting Functions   //////////////////////////////////////////


bool	SaveRaster2ImageFile(Raster<uchar>&image, const char *pname_file, const int dpi = 300);
int		LoadImageFile2Raster(Raster<uchar>&image, const char *pname_file);

#endif
